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Abstract 

In this paper, we consider the optimal design of photonic crystal band structures for two- 
dimensional square lattices. The mathematical formulation of the band gap optimization prob- 
lem leads to an infinite-dimensional Hermitian eigenvalue optimization problem parametrized by 
the dielectric material and the wave vector. To make the problem tractable, the original eigen- 
value problem is discretized using the finite element method into a series of finite-dimensional 
eigenvalue problems for multiple values of the wave vector parameter. The resulting optimization 
problem is large-scale and non-convex, with low regularity and non-differentiable objective. By 
restricting to appropriate eigenspaces, we reduce the large-scale non-convex optimization prob- 
lem via reparametrization to a sequence of small-scale convex semidefinite programs (SDPs) 
for which modern SDP solvers can be efficiently applied. Numerical results are presented for 
both transverse magnetic (TM) and transverse electric (TE) polarizations at several frequency 
bands. The optimized structures exhibit patterns which go far beyond typical physical intuition 
on periodic media design. 

1 Introduction 

The propagation of waves in periodic media has attracted considerable interest in recent years. 
This interest stems from the possibility of creating periodic structures that exhibit band gaps 
in their spectrum, i.e., frequency regions in which the wave propagation is prohibited. Band 
gaps occur in many wave propagation phenomena including electromagnetic, acoustic and elastic 
waves. Periodic structures exhibiting electromagnetic wave band gaps, or photonic crystals, have 
proven very important as device components for integrated optics including frequency filters 
waveguides [lOj, switches [20], and optical buffers [27j. 

The optimal conditions for the appearance of gaps were first studied for one-dimensional crys- 
tals by Lord Rayleigh in 1887 [18j. In a one-dimensional periodic structure, one can widen the 
band gap by increasing the contrast in the refractive index and difference in width between the 
materials. Furthermore, it is possible to create band gaps for any particular frequency by changing 

*This research has been supported through AFOSR grant FA9550-08- 1-0350 and the Singapore-MIT Alliance. 
^National University of Singapore, Center for Singapore-MIT Alliance, Singapore 117576, email: men@nus . edu . sg 
■^MIT Department of Aeronautics and Astronautics, 77 Massachusetts Ave., Cambridge, MA 02139, USA, email: 
cuongngQmit . edu 

§MIT Sloan School of Management, 50 Memorial Drive, Cambridge, MA 02142, USA, email: rfreund@mit.edu 
^MIT Department of Electrical Engineering and Computer Science, 77 Massachusetts Ave., Cambridge, MA 02139, 

USA, email: parriloQmit . edu 

"MIT Department of Aeronautics and Astronautics, 77 Massachusetts Ave., Cambridge, MA 02139, USA, email: 

peraire^mit . edu 



the periodicity length of the crystal. Unfortunately, however, in two or three dimensions one can 
only suggest rules of thumb for the existence of a band gap in a periodic structure, since no rigorous 
criteria have yet been determined. This made the design of two- or three-dimensional crystals a trial 
and error process, being far from optimal. Indeed, the possibility of two- and three-dimensionally 
periodic crystals with corresponding two- and three-dimensional band gaps was not suggested until 
100 years after Rayleigh's discovery of photonic band gap in one dimension, by Yablonovitch [25] 
and John [14j in 1987. 

From a mathematical viewpoint, the calculation of the band gap reduces to the solution of an 
infinite-dimensional Hermitian eigenvalue problem which is parametrized by the dielectric function 
and the wave vector. In the design setting, however, one wishes to know the answer to the question: 
which periodic structures, composed of arbitrary arrangements of two or more different materials, 
produce the largest band gaps around a certain frequency? This question can be rigorously ad- 
dressed by formulating an optimization problem for the parameters that represent the material 
properties and geometry of the periodic structure. The resulting problem is infinite-dimensional 
with an infinite number of constraints. After appropriate discretization in space and consideration 
of a finite set of wave vectors, one obtains a large-scale finite-dimensional eigenvalue problem which 
is non-convex and is known to be non-differentiable when eigenvalue multiplicities exist. The cur- 
rent state-of-the-art work done on this problem falls into two broad categories. The first kind tries 
to find the "optimal" band structure by parameter studies - based on prescribed inclusion shapes 
(e.g., circular or hexagonal inclusions) [9J or fixed topology [26j. The second kind attempts to use 
formal topology optimization techniques |T9l [71 [H [15] . Both approaches typically use gradient- 
based optimization methods. While these methods are attractive and have been quite successful 
in practice, the optimization processes employed explicitly compute the sensitivities of eigenval- 
ues with respect to the dielectric function, which are local subgradients for such non-differentiable 
problem. As a result, gradient-based solution methods often suffer from the lack of regularity of 
the underlying problem when eigenvalue multiplicities are present, as they typically are at or near 
the solution. 

In this paper we propose a new approach based on semidefinite programming (SDP) and sub- 
space methods for the optimal design of photonic band structure. In the last two decades, SDP has 
emerged as the most important class of models in convex optimization; see [TJ [21 [161 [221 [2l] . SDP 
encompasses a huge array of convex problems as special cases, and is computationally tractable 
(usually comparable to least-square problems of comparable dimensions). There are three distinct 
properties that make SDP very suitable for the band gap optimization problem. First, the un- 
derlying differential operator is Hermitian and positive semidefinite. Second, the objective and 
associated constraints involve bounds on eigenvalues of matrices. And third, as explained below, 
we can approximate the original non-convex optimization problem by a semidefinite program for 
which SDP can be well applied, thanks to its efficiency and robustness of handling this type of 
spectral objective and constraints. 

In our approach, we first reformulate the original problem of maximizing the band gap between 
two consecutive eigenvalues as an optimization problem in which we optimize the gap in eigenvalues 
between two orthogonal subspaces. The first eigenspace consists of eigenfunctions corresponding 
to eigenvalues below the band gap, whereas the second eigenspace consists of eigenfunctions whose 
eigenvalues are above the band gap. In this way, the eigenvalues are no longer present in our 
formulation; however, like the original problem, the exactly reformulated optimization problem 
is large-scale. To reduce the problem size, we truncate the high-dimensional subspaces to only 
a few eigenfunctions below and above the band gap |5l |T7], thereby obtaining a new small-scale 
yet non-convex optimization problem. Finally, we keep the subspaces fixed at a given decision 
parameter vector and use a reparametrization of the decision variables to obtain a convex semidef- 
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inite optimization problem for which SDP solution methods can be effectively applied. We apply 
this approach to optimize band gaps in two-dimensional photonic crystals for either the transverse 
magnetic (TM) or the transverse electric (TE) polarizations. 

The rest of the paper is organized as follows. In Section [2] we introduce the governing differential 
equations and the mathematical formulation of the band gap optimization problem. We then discuss 
the discretization process and present the subspace restriction approach. In Section [3] we introduce 
the semidefinite programming formulation of the band structure optimization, and lay out the 
optimization steps involved in solving the problem. Numerical results are presented in Section [4] 
for both the TE and TM polarizations in square lattices. Finally, in Section [5] we conclude with 
several remarks on anticipated future research directions. 



2 The Band Gap Optimization Problem 
2.1 Governing Equations 

Our primary concern is the propagation of electromagnetic linear waves in periodic media, and 
the design of such periodic structures, or photonic crystals, to create optimal band gaps in their 
spectrum. The propagation of electromagnetic waves in photonic crystals is governed by Maxwell's 
equations. The solutions to these equations are in general very complex functions of space and 
time. Due to linearity however, it is possible to separate the time dependence from the spatial 
dependence by expanding the solution in terms of harmonic modes - any time- varying solution can 
always be reconstructed by a linear combination of these harmonic modes using Fourier analysis. 
By considering only harmonic solutions, the problem is considerably simplified since it reduces to 
a series of eigenvalue problems for the spatially varying part of the solutions (eigenfunctions) and 
the corresponding frequencies (eigenvalues). 

In the absence of sources and assuming a monochromatic wave, i.e., with magnetic field (r, t) = 
H{r)e-'^\ and electric field E{r,t) = E{r)e-'''\ Maxweh's equations can be written in the fol- 
lowing form: 

X (V X E(r)) = (-\^E(r), in R^ 
e[r) \cJ 

where c is the speed of light, and e{r) is the dielectric function. In two dimensions, there are two 
possible polarizations of the magnetic and electric fields. In TE (transverse electric) polarization, 
the electric field is confined to the plane of wave propagation and the magnetic field H — (0, 0, i7) 
is perpendicular to this plane. In contrast, in TM (transverse magnetic) polarization, the magnetic 
field is confined to the plane of wave propagation and the electric field E — (0, 0, E) is perpendicular 
to this plane. In such cases, the Maxwell's equations can be reduced to scalar eigenvalue problems 

TE: -V- (^Vif(r)) = (5^)'if(r), in M^, (1) 

TM: -V ■{VE{r))= e{r)E{r), in . (2) 

Note that the reciprocal of the dielectric function is present in the differential operator for the TE 
case, whereas the dielectric function is present in the right-hand side for the TM case. 
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Figure 1: Left: A photonic crystal on a square lattice. The dashed box represents the primitive 
unit cell (fi), where a is the periodicity length of the lattice. Right: The reciprocal lattice, and 
the dashed box represents the first Brillouin zone (B). The irreducible zone is the green triangular 
wedge, and its boundary is denoted by dB. 
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For two-dimensional square lattices the dielectric function satisfies e{r) = £{r + JR), where R 
are the crystal lattice vector^ By applying the Bloch-Floquet theory [HIE] for periodic eigenvalue 
problems we obtain that 



H{r)^e'''-^Hk{r), and E{r) ^ e'''-^ Ek{r), 
where Hk{r) and Ej^{r) satisfy 



TE: (V + zfc)- ( ^{V + ik)Hk{r)^ = (^^y Hk{r), in f^, (3) 



TM : (V + ik) • ((V + ik)Ek{r)) = (^^^ s{r)Ek{r), in f^, (4) 

respectively. Thus, the effect of considering periodicity is reduced to replacing the indefinite periodic 
domain by the unit cell and V by V + ifc in the original equation, where fc is a wave vector 
in the first Brillouin zone B. Note that the unit cell and the Brillouin zone B depend on the 
lattice type (e.g., square or triangular lattices) as well as the crystal lattice vectors R. If we further 
take into consideration the symmetry group of the square lattice [23j, we only need to consider all 
possible wavevectors k on the irreducible Brillouin zone, or (under certain conditions) its boundary 
|13| . Figure [l] shows an example of the unit cell and the Brillouin zone for a square lattice. 
For notational convenience, we write the above equations in the following operator form 

Au = \Mu, in fi, (5) 

where, for the TE case, u = H^^r)^ A = u^^/c^^ and 

A{£, k) = -(V + ik) • (^(V + ik)^ , M = I; (6) 

whereas, for the TM case, u = Ek{r)^ A = u^y^/c^^ and 

A{k) = -(V + ik) • (V + ik), M{e) = e{r)L (7) 

Here / denotes the identity operator. We denote by (i^^, A^) the m-th pair of eigenfunction and 
eigenvalue of ([s]) and assume that these eigenpairs are numbered in ascending order: < A^ < 
A2 < ... < A^. 



^For a square lattice, R denotes the vectors spanned by {aei,ae2}, where ei and 62 are the unit basis vectors 
and a is the periodicity length of the crystal [13^ . 
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2.2 The Optimization Problem 

The objective in photonic crystal design is to maximize the band gap between two consecutive 
frequency modes. Due to the lack of fundamental length scale in Maxwell's equations, it can be 
shown that the magnitude of the band gap scales by a factor of s when the crystal is expanded by 
a factor of Therefore, it is more meaningful to maximize the gap-midgap ratio instead of the 
absolute band gap [13j. The gap-midgap ratio between and A^+^ is defined as 

_ i^hedB A^+^(g(r), fc) - sup^^^^ A^(g(r), fc) 

where dB represents the irreducible Brillouin zone boundary; see Figure [l] for example. 

A typical characterization of the dielectric function £{r) is the distribution of two different 
materials. Suppose that we are given two distinct materials with dielectric constants ^min and 5max 
where ^min < ^max- We wish to find arrangements of the materials within the unit cell Q which 
result in maximal gap-midgap ratio. To this end, we decompose the unit cell into disjoint 
subcells Ki^l < i < N^, such that ft = U^-^Ki and Ki n = for i j. Here we take this 
subcell grid to be the same as the finite element triangulation of the unit cell as we are going to 
discretize the continuous eigenvalue problem by the finite element method. Our dielectric function 
£{r) takes a unique value between ^min and ^max on each subcell, namely, £(r) = G M on 
and £min ^ £i ^ ^max- Howcvcr, duc to the symmetry of square lattice, we only need to define 
the dielectric function £{r) over part of the unit cell (1/8 of the unit cell). Hence, in general, the 
dielectric function e{r) is discretized into a finite dimensional vector s = (^1, • • • , ^n^) ^ (with 

^ ^e) which resides in the following admissible region: 

Qad = = (^1, • • • , ^ne) ^ IK""' • ^min < £i < ^max, 1 < ^ < ^s}- 

This region consists of piecewise-constant functions whose value on every subcell varies between ^min 
and £:max- Moreover, to render this problem computationally tractable, we replace the irreducible 
Brillouin zone boundary dB by a finite subset 

Sn, ^{ktedB, l<t< n J, 

where kt,l < t < Uk, are wave vectors chosen along the irreducible Brillouin zone boundary. As a 
result, the band gap optimization problem that maximizes the gap-midgap ratio between A^ and 
A^^-^ can be stated as follows: 

X ^^^keSn, A^+i(s, k) - maxfce^., A^(s, k) 

max J (s) = —, ^4 — — TT — — ^4 — - — — 

mmj^^Sr, A^+i(s, k) + max^^^., A^(s, k) 

(8) 

s.t. A(s^k)u^ = M{e)u^ ^ j = m,m + l, k ^ Sn^^ 

^min ^ ^ ^max? 1 ^ ^ ^ '^e- 

In this problem a subtle difference between TE and TM polarizations lies in the operators of the 
eigenvalue problem: A and Ai take the form of either ([g]) for the TE case or for the TM case. 
In either case, note that the eigenvalue problems embedded in ([s]) must be addressed as part of any 
computational strategy for the overall solution of Q. 

2.3 Discretization of the Eigenvalue Problem 

We consider here the finite element method to discretize the continuous eigenvalue problem ([5]). 
This produces the following discrete eigenvalue problem 

Ah{e, k)ul = XlMh{e)ul j = 1, . . . , AA, G 5,,, (9) 
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where k) G C-^^-^ is a Hermitian stiffness matrix and M/^(s) G R-^^-^ is a symmetric positive 

definite mass matrix. These matrices are sparse and typicahy very large (AA ^ 1). We consider 
the approximate eigenvalues in ascending order: < < • • • < A"^. 

It is important to note that the dependence of the above matrices on the design parameter 
vector e is different for the TE and TM polarizations. In the TE case, A™ depends on e and 
does not, whereas in the TM case depends on e and ^™ does not. More specifically, since 

e{r) is a piecewise-constant function on fi, the s-dependent matrices can be expressed as 

A^is, fc) = f; Ul^ik), M^is) = |:..M™, (10) 

i=i ^ i=i 

where the matrices A^^{k) and M™, I < i < are independent of e. We note that A^^{e, k) 
is linear with respect to 1/^^, 1 < z < n^, while M™(s) is linear with respect to 6^, 1 < z < 



He. The affine expansion (10) is a direct consequence of the fact that we use piecewise-constant 
approximation for the dielectric function £{r). (In the TE case, we will shortly change our decision 
variables to = 1/^^, 1 < z < n^, so as to render A™ affine in the variables yi, . . . , yue-) 

After discretizing the eigenvalue problem ([s]) by the finite element method, we obtain the 
following band gap optimization problem: 

UCidbX. J J^yS J — ^pj 

maxfce^., Xh ^) + ^^^keSn^ ^hi^^ ^) 

(11) 

s.t. Ah{e,k)ul = XlMh{e)ul, j = m,m + l, k e Snj,, 

^min ^ ^ ^max? 1 ^ ^ ^e- 

Unfortunately, this optimization problem is non-convex; furthermore it suffers from lack of regu- 
larity at the optimum. The reason for this is that the eigenvalues A^ and A^^^ are typically not 
smooth functions of e at points of multiplicity, and multiple eigenvalues at the optimum are typical 
of structures with symmetry. As a consequence, the gradient of the objective function J(s) with 
respect to e is not well-defined at points of eigenvalue multiplicity, and thus gradient-based descent 
methods often run into serious numerical difficulties and convergence problems. 



3 Band Structure Optimization 

In this section we describe our approach to solve the band gap optimization problem based on a 
subspace method and semidefinite programming (SDP). In our approach, we first reformulate the 
original problem as an optimization problem in which we aim to maximize the band gap obtained by 
restriction of the operator to two orthogonal subspaces. The first subspace consists of eigenfunctions 
associated to eigenvalues below the band gap, and the second subspace consists of eigenfunctions 
whose eigenvalues are above the band gap. In this way, the eigenvalues are no longer explicitly 
present in the formulation, and eigenvalue multiplicity no longer leads to lack of regularity. The 
reformulated optimization problem is exact but non-convex and large-scale. To reduce the problem 
size, we truncate the high-dimensional subspaces to only a few eigenfunctions below and above the 
band gap [S^^TT], thereby obtaining a new small-scale yet non-convex optimization problem. Finally, 
we keep the subspaces fixed at a given decision parameter vector to obtain a convex semidefinite 
optimization problem for which SDP solution methods can be efficiently applied. 
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3.1 Reformulation of the Band Gap Optimization Problem using Subspaces 

We first define two additional decision variables: 

mm A^+^e, k) , Af max A^(£, k) , 



and then rewrite the original problem (11) as 



Pq : max 
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(12) 



s.t. \^{e, k) < \i , XI < A^+H^, k), Vfc e 5„ 

fcX = A;^M;,(£)n}^, Vfc e Sn 

Ah{e, k)u^+' = A-+iM;,(£)n-+\ Vfc e 5„ 

^min ^ ^ "S^max? ^ — • • • 5 

K . > 0. 

Next, we introduce the following matrices: 

$-(fc) := Mik) I Kik)] ■■= [uUe,k) ...u^{e,k) \ u^+\s,k) . . . uff {e, k)], 

where $f (fc) and ^u(k) consist of the first m eigenvectors and the remaining J\f — m eigenvectors, 
respectively, of the eigenvalue problem: 

Ah{e, k)ul = XlMh{e)ul l<j<N. 

We will also denote the subspaces spanned by the eigenvectors of $f (fc) and as sp($|(fc)) 

and sp($^(fc)), respectively. 

The first three sets of constraints in ( [T2| ) can be represented exactly as 



<^r{k)[Ahie, k) - XiMhieWdk) ^0, Vfc G 5, 
^l*{k)[Ah{s,k) - XlMhisWuik) hO, Vfc G 5, 



where is the Lowner partial ordering on symmetric matrices, i.e., A^Bif and only if A — B 
is positive semidefinite. We therefore obtain the following equivalent optimization problem: 

Pi : max — j- 

^MM K + K 

s.t. ^*{k)[Ah{B,k)-XiMh{s)]mk)<{), ykeSn,, (13) 
^l*{k)[Ahie,k) - XlMhie)]^l{k) hO, Vfc e Sn„ 

^min ^ ^ "S^max? ^ — • • • 5 

K . > 0. 

Although the reformulation Pi is exact, there is however a subtle difference in the interpretation 
of Pq and Pi: Pq can be viewed as maximizing the gap-midgap ratio between the two eigenvalues 
and X^^^; whereas Pi can be viewed as maximizing the gap-midgap ratio between the two 
subspaces sp($|(fc)) and sp($!^(fc)). The latter viewpoint allows us to develop an efficient subspace 
approximation method for solving the band gap optimization problem as discussed below. 
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3.2 Subspace Approximation and Reduction 

Let us assume that we are given a parameter vector e. We then introduce the associated matrices 

Mik) I ^lik)] = [ul{e,k) ...u^{e,k) \ u^+\e,k) ...u^{e,k)] , 

where $f (A;) and $^(A;) consist of the first m eigenvectors and the remaining M — m eigenvectors, 
respectively, of the eigenvalue problem 

Ahie, k)ui = X{Mh{e)ul l<j<N. 

Under the presumption that sp($|(fe)) and sp(#^(fc)) are reasonable approximations of sp(#|(fc)) 
and sp($^(fc)) for s near e, we replace $f(fc) with and with to obtain 

Pi : max ~ 



s.t. $f (fc)[A(£,fc)-AlM^(£)]#f(fe) ^0, VfeG<Sn,, (14) 

^min ^ ^2 ^ ^max? ^ ~ • • • ? 

K > > 0. 

Note in ^^at the subspaces sp($|(fc)) and sp($^(fc)) are approximations of the subspaces 
sp($|(fc)) and sp($!^(fc)) and are no longer functions of the decision variable vector e. 

Note also that the semidefinite inclusions in P2 are large-scale, i.e., the rank of either the first 
or second inclusion is at least A/'/2, for each k G 5^^;., and J\f will typically be quite large. In order 
to reduce the size of the inclusions, we reduce the dimensions of the subspaces by considering only 
the "important" eigenvectors among u\{e^k) ... fc), fc) ... namely those 

ttfc eigenvectors whose eigenvalues lie below but nearest to AJ[^(s,fc) and those hy. eigenvectors 
whose eigenvalues lie above but nearest to A^^^(s, fc), for small values of a^, 6^, typically chosen 
in the range between 2 and 5, for each k G 5^^. This yields reduced matrices 

n,+H(k) := m^{k) I = ...u^{s,k) I u^+\s,k) ... 

Substituting ^a^i^) place of and ^f^(fc) in place of in the formulation P| yields 

the following reduced optimization formulation: 

\u _ \i 

Pi : max ^ 

S.t. $f*(fc)[A(s,fc) - XiMnieMAk) ^0, Vfc G 5,,, (15) 
$g(fc)[A,(s,fc) - XlMH{e)]^^{k) ^0, Vfc G 5,,, 

^min ^ ^ ^max? ^ ~ 1? • • • 5 

a;: , Al > 0. 

In this way the formulation seeks to model only the anticipated "active" eigenvalue con- 
straints, in exact extension of active-set methods in nonlinear optimization. The integers a^, are 
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determined indirectly through user-defined parameters ri > 0, and Tu > 0, where we retain only 
those eigenvectors whose eigenvalues are within lOOr/% beneath X^(s^ k) or whose eigenvalues are 
within 100r^i% above X^^^{e,k). This translates to choosing a^^bk G N+ as the smallest integers 
that satisfy 

A^(s,fc)-Ar"'°+\^,fc) < ^ Ang,fc)-Ar°'°(g,fc) 

\^{e,k) - - A-(£,fc) 

K^\£,k) - - X^+\s,k) 

The dimensions of the resulting subspaces sp($afc(fc)) and sp($^^(fc)) are typically very small 
6fc ^ 2, . . . , 5). Furthermore, the subspaces are well-spanned by including all relevant eigenvec- 
tors corresponding to those eigenvalues with multiplicity at or near the current min/max values. 

We observe that has significantly smaller semidefinite inclusions than if the full subspaces 
were used. Also, the subspaces are kept fixed at e in order to reduce the nonlinearity of the 
underlying problem. Furthermore, we show below that for the TE and TM polarizations that P3 
can be easily re-formulated as a linear fractional semidefinite program, and hence is solvable using 
modern interior-point methods. 

3.3 Fractional SDP Formulations for TE and TM Polarizations 

We now show that by a simple change of variables for each of the TE and TM polarizations, problem 
P3 can be converted to a linear fractional semidefinite program and hence can be further converted 
to a linear semidefinite program. 

TE polarization. We introduce the following new decision variable notation for convenience: 

y := (l/i,I/2,...,I/nJ := (l/£:i,...,l/£ne,A^,A)J) , 
and set ^min = 1/^max and ^max = 1/^min- We also amend our notation to write various functional 



dependencies on y instead of e such as etc. Utilizing (10), we re- write for the TE 

polarization as 



Prpg : max 



Vny HUy-l 



y Vny + Vny-l 



s.t. 



EZf y^AEik) - yny-iM^^ < {k) ^0, Vfc G 5. 



rife 5 



(16) 



Uinm ^ Vi ^ I/max 5 



TE 



i = 1, . . . , 



VUy — 1 5 VUy 



> 0. 



We note that the objective function is a linear fractional expression and the constraint functions 
are linear functions of the variables y. Therefore is a linear fractional SDP. Using a standard 
homogenization a linear fractional SDP can be converted to a linear SDP0 

^Indeed, for notational simplicity consider a linear fractional optimization problem of the form max^ ^rf subject 
to b — Ax G Ki, X G K2, where d^x > for all feasible x and Ki, K2 are convex cones. Then this problem is 
equivalent to the problem max^^,6' c^'^ subject to bO — Aw G i^i, it; G i^2, d^w = 1, ^ > 0, under the elementary 
transformations x ^ (w/0) and [w^O) ^ (x/d^x, see [6l[8]. 
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TM polarization. We introduce slightly different decision variable notation for convenience: 

z := (2:1,2:2, ... ,2:^J := (^i, . . . , , l/A^, 1/A)J), 

and set Zmin = ^min and 2:max = ^max- Similar to the TE case, we amend our notation to write 
various functional dependencies on z instead of e such as $f (fc), etc. Noting that 



-1 + 



utilizing (10), and multiplying the semidefinite inclusions of (15) by 2:^^-1 and 2:^^, respectively, we 



re-write P3 for the TM polarization as 



TM 



max 

z 



s.t. 



^min ^ Zi ^ -2:max 



i 



1, ■ . . ,n. 



(17) 



Here again the objective function is a linear fractional form and the constraint functions are linear 
functions of the variables z. Therefore ^ linear fractional SDP with format similar to that 

Since both P^^ and Pr^y[ are linear fractional semidefinite programs, they can be solved very 
efficiently by using modern interior point methods. Here we use the SDPT3 software |21| for this 
task. 



3.4 Main Algorithm 

We summarize our numerical approach for solving the band gap optimization problem of the TE 
polarization in the following table. Essentially the same algorithm (with the modifications described 
in the previous section) is used to solve the band gap optimization problem of the TM polarization. 



4 Results and Discussions 
4.1 Model Setup 

We consider a two-dimensional photonic crystal confined in the computational domain of a unit cell 
of the square lattice, and with square domain = [— 1, 1] x [— 1, 1]. The domain is decomposed 
into a uniform quadrilateral (in particular, we use square elements for the square lattice) grid of 
dimensions 64 x 64, which yields a mesh size of /i = 1/32 and 4096 linear square elements. 

The dielectric function e is composed of two materials with dielectric constants ^min = 1 (air) 



and Smax = 11-4 (GaAs). As mentioned earlier in Section 2.2, the symmetry of the lattice can be 
exploited to further reduce the dielectric function to be defined in only 1/8 of the computation 
domain. The number of decision variables relating the dielectric material (s^, i = 1, 2, . . . ^rie) is 
thus reduced to = (1 + 32) x 32/2 = 528. Figure [2] shows an illustration of a coarse mesh 
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Implementation Steps 


Step 


1. 


Start with an initial guess and an error tolerance €tob s-nd set y := y^. 


Step 


2. 


For each wave vector k G Sn,., do: 






Determine the subspace dimensions and 6^. 






Compute the matrices $afe(fc) and 


Step 


3. 


Form the semidefinite program P^^- 


Step 


4. 


Solve for an optimal solution y*. 


Step 


5. 


U \\y* — y\\ < etoi, stop and return the optimal solution y* . 






Else update y ^ y* and go to Step 2. 



Table 1: Main algorithm for solving the band gap optimization problem. 



(16 X 16) and dielectric function for the square lattice to aid visualization; note that the actual 
computational mesh (64 x 64) is finer than this one. The shaded cells represent those modeled by 
6:, and the rest are obtained through symmetry. Furthermore, in this case, the irreducible Brillouin 
zone B is the triangle shown in Figure [TJ with rtk — 12 fc-points taken along the boundary of this 
region {dB). Band diagrams plotted in the figures below show the eigenvalues moving along the 
boundary of 0, from F to X to M and back to F. 



x:::n::n 



Figure 2: An illustration of a coarse mesh (16 x 16) and dielectric function for the square lattice. 
The shaded cells indicate the decision variables relating the dielectric material (5^, z = 1, 2, . . . , n^). 
Note that the actual computational mesh (64 x 64) is finer than this one. 



4.2 Choices of Parameters 
4.2.1 Initial configuration 

Because the underlying optimization problem may have many local optima, the performance of 
our method can be sensitive to the choice of the initial values of the decision variables y^, which 
in turn depend on the initial configuration . Indeed, different initial configurations do lead to 
different local optima as shown in Figure |3] for the second TE band gap and in Figure |4] for the 
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fourth TM band gap. Therefore, the choice of the initial configuration is important. We examine 
here two different types of initial configurations: photonic crystals exhibiting band gaps at the low 
frequency spectrum and random distribution. 

The well-known photonic crystals (e.g., dielectric rods in air - Figure |4]^a), air holes in dielectric 
material, orthogonal dielectric veins - Figure |3]^d)) exhibit band-gap structures at the low frequency 
spectrum. Such a distribution seems to be a sensible choice for the initial configuration as it 
resembles various known optimal structures [4J. When these well-known photonic crystals are 
used as the initial configuration, our method easily produces the band-gap structures at the low 
frequency mode (typically, the first three TE and TM modes). On the other hand, maximizing the 
band gap at the high frequency mode (typically, above the first three TE and TM modes) tends 
to produce more complicated structures which are very different from the known photonic crystals 
mentioned above. As a result, when these photonic crystals are used as the initial configurations 
for maximizing the band gap at the high frequency mode, the obtained results are less satisfactory. 

Random initial configurations such as Figures lll^a) and|4];d)) have very high spatial variation 
and may thus be suitable for maximizing the band gap at the high frequency mode. Indeed, we 
observe that random distributions often yield larger band gaps (better results) than the known 
photonic crystals for the high frequency modes. Of course, the random initialization does not 
eliminate the possibility of multiple local optima intrinsic to the physical problem. In view of this 
effect, we use multiple random distributions to initialize our method. In particular, we start our 
main algorithm with a number of uniformly random distributions as initial configurations to obtain 
the optimal structures in our numerical results discussed below. 




(a) Initial crystal configuration #1 (b) Optimized crystal structure #1 (c) Optimized band structure #1 




(d) Initial crystal configuration #2 (e) Optimized crystal structure #2 (f) Optimized band structure #2 
Figure 3: Two locally optimal band gaps between and in the square lattice 
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(a) Initial crystal configuration #1 (b) Optimized crystal structure #1 (c) Optimized band structure #1 




(d) Initial crystal configuration #2 (e) Optimized crystal structure #2 (f) Optimized band structure #2 
Figure 4: Two locally optimal band gaps between A^j^ and X^y^ in the square lattice 

4.2.2 Subspace dimensions 

The dimensions of the subspaces sp($afc(fc)) and sp($^^(fc)) are determined indirectly by the 
parameters n and Tu- A good choice of n (and ru) is one that returns <^ Af (and <^ A/"), and 
at the same time includes the "important" eigenvectors to enhance convergence to an optimum. In 
our numerical experiments, we choose r^^ = r/ = 0.1 which in turn leads to the resulting subspace 
dimensions and in the range of [2, 5]. Moreover, we find that choosing larger values of Tu and 
ri (e.g., ri — Tu — 0.2), which in turn increases and b^ and hence increases computational cost, 
does not yield fewer iterations than choosing r-^ = r/ = 0.1. 

4.3 Computational Cost 

With all the programs implemented in MATLAB and the computation performed on a Linux PC 
with Dual Core AMD Opteron 270, 1.99GHz, a successful run of the algorithm can typically be 
done in 2-30 minutes including 5-30 outer iterations, i.e., passes of Steps 2-5 of the main algorithm 
in Table [l} An example of the computational cost and outer iterations for different band gap 
optimization is shown in Table [2] as a general illustration of our computational experience. 

We point out that these numbers merely represent one set of possibilities; variations in the 
numerical results are likely to occur with different random initial configurations. Nevertheless, the 
computation cost does serve as an indication of the general level of difficulty of finding a solution 
in each problem. In general, lower eigenvalue band gap optimization problems are easier to solve 
(at least to local optima). Moreover, the table illustrates that TM problems usually solve faster 
and require fewer outer iterations. This latter observation is consistent with the result reported in 
[15], and is possibly due to the high non-convexity of the original TE optimization problem. 

Before ending this section, we point out some possible ways to improve the computational cost 
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AA^;^' 


AAII' 


AA'^:t' 


AA^:t' 


AA'ff 




^^9,10 


^^10,11 


Execution time (min) 


5.7 


2.5 


8.9 


20.4 


17.9 


20.5 


19.4 


27.3 


26.4 


25.8 


Outer Iterations 


11 


8 


29 


26 


18 


25 


15 


27 


19 


23 






AA^:^ 


AA^:f 


AAi:r^ 


AA'^:^^ 


AA^;^r^ 


AA'^;r^ 




^^9,10 


^^10,11 


Execution time (min) 


1.8 


5.6 


3.5 


5.4 


11.7 


9.5 


10.8 


3.9 


11.2 


9.5 


Outer Iterations 


4 


9 


5 


7 


16 


9 


9 


9 


12 


10 



Table 2: Example of computation time and the number of outer iterations of a successful run 
for optimizing various band gaps, for both TE and TM polarization. Here AXf^-^ denotes the 
gap-midgap ratio between the i^^ and (z + 1)^* eigenvalue for the TE polarization. 





(a) (b) 

Figure 5: (a) Random starting structure with translation, rotation, and reflection symmetry, 3x3 
unit cells in square lattice, (b) Band structure before optimization. 



of our procedure. For the eigenvalue calculation, it is probably helpful to apply a more efficient 
eigensolver (we used MATLAB's eigs function in the current implementation). Another promising 
approach is to explore mesh adaptivity and incorporate non-uniform meshing for the representation 
of the dielectric function, as well as the eigenvalue calculation. As further discussed in Section [5j 
it should be possible to significantly reduce the number of decision variables and computation cost 
with this approach. 



4.4 Optimal Structures 

We start the optimization procedure with a random distribution of the dielectric, such as the one 
shown in Figure [5]^ a) . The corresponding band structures of the TE and TM fields are shown 
in Figure [5]^b). In Figure [6j we present an example of the evolution of the crystal structure as 
the optimization process progresses. (The light color indicates the low dielectric constant and the 
dark color denotes the high dielectric constant.) As illustrated in Figure [7| the gap-midgap ratio 
starts from a negative value (—8.93%) corresponding to the random configuration (Figure [t]^ a)) and 
increases up to +43.90% corresponding to the optimal configuration (Figure [7]^f)) at which time 
the optimization process terminates successfully. Another example of the optimization evolution 
for TE polarization is shown in Figure [8] and Figure |9j in which the gap-midgap ratio increases 
from -39.21% to +29.23%. 
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(a) 



(b) 



(c) 



(d) 



(e) 



(f) 



Figure 6: The evolution of the square lattice crystal structure for optimizing the gap-midgap ratio 
between A^|^ and A^j^. 






(a) 



(b) 




X % = 32.67% 




(d) 



(e) 



(f) 



Figure 7: The corresponding band structure (of Figure [g])) and the gap-midgap ratio between A^j^ 
and A|.j^ in the square lattice. 



15 





16 



(a) Optimal crystal structure 



(b) Optimal band structure 



Figure 10: Optimization of band gap between A^|^ and A^|^ in the square lattice. 




(a) Optimal crystal structure (b) Optimal band structure 

Figure 11: Optimization of band gap between A^|^ and A^|^ in the square lattice. 



In Figures [10] through [T9| we present only plots of the final optimized crystal structures and 
the corresponding band structures for the 6*^ through 10*^ optimized band gaps for TE and TM 
polarizations. We see that the optimized TM band gaps are exhibited in isolated high-5 structures, 
while the optimized TE band gaps appear in connected high-£: structures. This observation has also 
been pointed out in [13j (p75) Hhe TM band gaps are favored in a lattice of isolated high-e regions, 
and TE band gaps are favored in a connected lattice^\ and observed in [15j previously. For both TE 
and TM polarizations, the crystal structures become more and more complicated as we progress 
to higher bands. It would be very difficult to create such structures using physical intuition alone. 
The largest gap-midgap ratio for the TM case is 43.9% between the seventh and eighth frequency 
bands, while the largest ratio for the TE case is 44.1%, also between the seventh and eighth bands. 
The results presented here are not guaranteed to be globally optimal, as pointed out in Section 
4.2.1. While most crystal structures in the TM cases appear similar to those presented in [15], we 
have shown quite different TE structures. A qualitative comparison between the two results in the 
background indicates larger band gaps (both in absolute value and in the gap-midgap ratio) in our 
results. 
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(a) Optimal crystal structure 



(b) Optimal band structure 



Figure 12: Optimization of band gap between A^^i and X^yv in the square lattice 




(a) Optimal crystal structure 



(b) Optimal band structure 



Figure 13: Optimization of band gap between A^|^ and A^^|^ in the square lattice 




(a) Optimal crystal structure 



(b) Optimal band structure 



Figure 14: Optimization of band gap between X}py^ and A^^|^ in the square lattice 
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(a) Optimal crystal structure (b) Optimal band structure 

Figure 15: Optimization of band gap between A^g and A^g in the square lattice. 
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(a) Optimal crystal structure (b) Optimal band structure 

Figure 18: Optimization of band gap between A^g and A^g in the square lattice. 
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5 Conclusions and Future Work 



We have introduced a novel approach, based on reduced eigenspaces and semidefinite programming, 
for the optimization of band gaps of two-dimensional photonic crystals on square lattices. Our 
numerical results convincingly show that the proposed method is very effective in producing a 
variety of structures with large band gaps at various frequency levels in the spectrum. 

Since our computational techniques make essential use of the finite element method, we antic- 
ipate that notions of mesh adaptivity can be easily incorporated into our approach, and thus its 
computational efficiency will be improved even further. For example, one can start with a relatively 
coarse mesh and converge to a near-optimal solution, and then judiciously refine the finite element 
mesh (e.g., refining elements at the interface of dielectric materials) using the current optimal solu- 
tion at the coarser mesh as the new initial configuration. We intend to explore this approach and 
report the details and results in a forthcoming paper. 

The main strengths of our proposed approach to solve eigenvalue gap optimization problem 
is the fact that SDP-based methods do not require explicit computation of (sub-)gradients of the 
objective function (which are ill-defined in the case of eigenvalue multiplicities), hence maintaining 
the regularity of the formulation. The approach proposed in this paper can also be readily extended 
to deal with more general problems, such as the optimization of photonic crystals in combined TE 
and TM fields, optimizing multiple band gaps, dealing with other types of lattices (e.g. triangular), 
as well as modeling and optimizing the design of three-dimensional photonic crystals. 
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